Variation of Inner Radius of Dust Torus in NGC4151 
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ABSTRACT 

The long-term optical and near infrared monitoring observations for a type 
1 active galactic nucleus NGC 4151 were carried out for six years from 2001 to 
2006 by using the MAGNUM telescope, and delayed response of flux variations 
in the K(2.2/j,m) band to those in the V(0.55/j,m) band was clearly detected. 
Based on cross correlation analysis, we precisely measured a lag time At for 
eight separate periods, and we found that At is not constant changing between 
30 and 70 d during the monitoring period. Since At is the light travel time 
from the central energy source out to the surrounding dust torus, this is the 
first convincing evidence that the inner radius of dust torus did change in an 
individual AGN. In order to relate such a change of At with a change of AGN 
luminosity L, we presented a method of taking an average of the observed U-band 
fluxes that corresponds to the measured value of At, and we found that the time- 
changing track of NGC 4151 in the At versus L diagram during the monitoring 
period deviates from the relation of At oc L 5 expected from dust reverberation. 
This result, combined with the elapsed time from period to period for which At 
was measured, indicates that the timescale of dust formation is about one year, 
which should be taken into account as a new constraint in future studies of dust 
evolution in AGNs. 

Subject headings: dust, extinction — galaxies: active — galaxies: individual (NGC 
4151) — galaxies: Seyfert — infrared: galaxies 
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Introduction 



The unified model of active galactic nuclei (AGNs) (jAntonuccilll993l ; lUrry and Padovali 



19951 ) assumes the existence of dust torus that surrounds the central hierarchical structure 
consisting of a super massive black hole, an accretion disk, and a broad line region 
(BLR). The dust grains in the torus absorb the UV/optical continuum emission from the 
accretion disk and re-radiate in the near infrared (IR) wavelength region with some lag 
time corresponding to the light travel time from the accretion disk to the inner radius 
of dust torus. Since the heated dust eventually sublimates at a constant temperature of 
1500-1800K, a more luminous AGN should have a larger dust torus and hence a larger 
l ag time, yielding a co rrelation of At oc L 5 between lag time At and AGN luminosity L 



(parvainis 



1987 



199J). 



In 


: act, 


based on the 


ong-t 


(Yoshii 


2002 




Yoshii et al. 


20031) 



ong-term multicolor monitoring data from the MAGNUM project 



Suganuma et al. 



(}2006|) recently 



and available archival data 
reported such a correlation from a sample of AGNs spanning a wide range of absolute 
^-magnitude from My = — 15 to —24. The optical luminosity is a good in dicator of UV 



lumin osity, because their variations are well synchronized with each other (jEdelson et al 



19961 ) . Therefore, the reported correlation in the V band indicates that the inner radius of 
dust torus scales with the UV luminosity as well, which provides a strong piece of evidence 
for dust reverberation proceeding in the central region of AGNs. 

Given the above result of scaling among AGNs of different luminosities, it is possible 
to examine how At would scale with changing L in an individual AGN. This is currently 
an only diagnosis on its inner structure which is too small to be resolved by direct imaging. 
However, since a range of flux variation in an AGN is as small as < 0.1 mag, detection of 
a change of At is indeed a challenge and needs an elaborate method that enables precise 
estimation of At. 
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We here attempt to detect, if any, a change of At for the case of a nearby Seyfert 1 
AGN of NGC 4151, using the six year log accurate and frequent monitoring data for bright 
AGNs in the program of MAGNUM. 

In Section [21 we show details of observation of NGC4151. In Section [3] and Section H] 
we formulate estimation of At and L, by making use of model light curves. In Section Owe 
apply the formulation to the case of NGC 4151 and derive its time-changing track in the At 
versus L diagram. In Section [6] we discuss time scale of dust formation in the central region 
of NGC 4151. 



Observations of NGC 4151 



NGC4151 is one of the bright AGNs most frequently observed since 2001 by the 
MAGNUM telescope. The monitoring period analyzed in this Letter covers a length of 
2000 d from MJD 52, 000 to 54, 000 (2001-2006), and photometric images were obtained for 
about 220 nights in the V and K ba nds. The observations , photometry, and data analysis 



are mostly in common with those in 



Minezaki et al. 



(120041 ). and we do not repeat them 



except for some important basics below. 

Each night we took images of the AGN and the reference star alternatively using the 
multicolor imaging photometer (MIP) with the viewing angle of 1'.5, and performed relative 
photometry of the AGN with an aperture size of (f) = 8". 3. This size was chosen to minimize 
the photometric error and the flux variation in the aperture by the seeing effect. 

The flux in the aperture contains the flux from the host galaxy, especially from its 
bulge. The host component was estimated by th e model fitting to its surface brightness as 



A/ gal = 1^.95 mJy (?) and f K g & \ = 44.22 mJy (IMinezaki et al. 



20041 ). and we subtracted 



these contributions from the V and K light curves, respectively. The narrow line flux 
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was also significant in the V band for NGC4151, so we subtracted 10.41 mJy following 
?. Furthermore, the IR flux contains the flux from the a ccretion disk, which makes At 
systematically shorter than the actual lag. According to 



Tomita et al 



(120061 ) . the IR 



disk component should vary to synchronize with the optical variation that is originated 
genuinely from the accretion disk. This component in the K band was then estimated as 
/W,acc = (yy/v K ) a f Vl and we also subtracted this contribution from the K light curve by 
assuming a = for the disk emission between the V and K bands. 

The resulting V and K light curves of NGC 4151 from 2001 to 2006 are shown in 
Figure[TJ T he V light curve shown he re is consistent with the optical continuum light curves 



reported by 



Shapovalova et al. 



(120081 ). We see that the optical flux of NGC 4151 reached 
a minimum around MJD 52, 000 and then increased until reaching a maximum at MJD 
52, 700. Thereafter, it turned to decrease again to a minimum at MJD 53, 500. Apparently, 



the flux variation in the K band followed that in the V band with clear time delay. 



3. Method of Analysis 

The light curve shown in Figure Q] is complex, with many peaks and valleys. From it we 
need to extract a table of At as a function of continuum flux which has sufficient precision 
for us to test the hypothesis that At is proportional to fy 5 as the flux varies, and/or to 
understand any difference from that expected result. 



The techniques use d previously are not adequate. In our series of papers on dust 



reverberation in AGNs ( Minezaki et al. 



2004 



Suganuma et al. 



2006 



Tomita et al. 



20061 ) 



we estimated At and its error based on t he Cross Correlation Fun ction (CCF) analysis by 



making use of the structure function (SF, I White &: Peterson 



1994 ) for the interpolation of 



optical flux variability during the monitoring period. In the equal sampling (ES) scheme of 



- 6- 



interpolation which we used for the CCF analysis (ISuganuma et al.ll2006l ). the number of 
data points generated at equal intervals of time can exceed that actually observed. In such 
a case, the simulated light curve shows larger flux variation compared to an overall trend 
seen from the observed light curve, especially when the data points are generated far from 
the location of observed epochs. This necessarily gives rise to too large an error of At. 

We here adopt the bi-directional interpolation (BI) scheme, where all l/-band flux data 
are paired with f^-band flux data and at least either V or K flux in each pair is measured 
by actual observation. While the BI scheme overcomes a problem of generating too many 
data points for CCF calculations, the CCF still shows an unreal peak with increasing 
the number of data pairs. This drawback is remedied by estimating At from the CCF 
centroid around the peak which is interprete d as a luminosity- weighted distance of the dust 



torus from the AGN central energy source (IKoratkar fc Gaskell 



19911 ). As a practice of 



determining the CCF range for centroid estimation, the threshold is usually set to be 0.8 
times the CCF-peak value. If the threshold is larger than 0.9, At is influenced by an unreal 
peak of the CCF. On the other hand, if it is smaller than 0.7, At is necessarily influenced 
by sub-features such as a flared foot of the CCF. 

We so far used a simple mean of observed V^-band fluxes as an indicator of optical 
luminosity that corresponds to At, but this could seriously be biased by the two factors. 
First, the flux mean depends on the sampling frequency of light curve during the monitoring 
period. When an AGN is often observed in the bright phase rather than in the faint phase, 
the simple flux mean preferentially gives a value of brighter flux. Second, the flux mean 
depends on the shape of light curve. When an AGN undergoes an abrupt burst on a flat 
light curve, the simple flux mean gives a value which is more or less the same as the value 
of flat part of the light curve, although the flat part contributes very little to the value of 
cross-correlation (CC) coefficient for the estimation of At. 
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Considering above, we take a mean of V-band fluxes generated at equal intervals of 
time, after excluding the data which statistically show no flux variation in time bins. To 
pick up such data, we use the equally interpolated optical light curve for which the size 
of time bins is taken as the mean interval of monitoring observations. The variability in 
each bin is tested by \ 2 with a significance level of 0.99, and the data which are judged 
as showing no flux variation are excluded in calculating the mean of V^-band fluxes. The 
distribution of this flux mean for all simulated light curves is then used to derive an average 
optical flux f v and its error. 

4. A Test Using Model Light Curves 

We next use model light curves to verify that our analysis procedure reliably finds a 
characteristic lag At and its associated flux fy for each of the individual variability events 
that combine to make the complex light curve shown in Figure 1. In order to quantify 
the above estimation of a characteristic lag At and its associated flux fy as well as their 
errors, we model the light curves to calculate CCFs. Given a V light curve fv(t), we 
construct a corresponding K light curve as fxif) = fv(t — At), where At/50 = (/y/200) 7 
for 7 = 0.5 ± 0.3 which roughly reproduce the lags seen in Figured! Here, the units are days 
in time and arbitrary in flux. Then we generate 30 data points along each of the V and K 
light curves, and disperse the points according to a Gaussian distribution of photometric 
errors taken from the actual observations by the MAGNUM telescope. 

The models include three characteristic shapes of (1) linear bottom, (2) quadratic 
bottom, (3) bursting Gaussian peak, as shown in order from top to bottom in the left 
column of Figure [2J The models 1 to 3 are shown in Figure [2] to examine the dependence 
on shape of time-symmetric incident V light curve used for CCF calculation. In addition to 
all these models with a standard index of 7 = 0.5, we consider two cases of 7 = 0.2 and 0.8 
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for each of light curve shapes 1 and 2. We note, however, that 7 > 0.5 is considered only 
for academic purpose and contradicts with the geometrical 1 / r 2 -dilution of the incident flux 
from the central energy source. 



Following 



Suganuma et al. 



(120061 ). we applied the BI method to each of model 
light curves to calculate the CCF and also the CC centroid distribution (histogram) by 
Monte-Carlo simulations, as shown in the middle column of Figure [2] for 7 = 0.5. The 
results of At and its error are given on the panels of this column and in Table [TJ The 
results of At for the models 1 to 3 give values around 80 d. This indicates that At could be 
estimated independently of the shape of light curve. 

Now we derive an average flux fy for each of the models 1 to 3, following the way 
described in Section [3J The horizontal arrows in the left column of Figure [2] specify 
the periods in which the data are rejected in the estimation of flux mean due to their 
invariability. We calculate the distribution of this flux mean for the light curves generated 
by Monte-Carlo simulations, and present it as a histogram from which fy and its error are 
derived, as shown in the right column of Figure [2] for 7 = 0.5. The results of fy and its 
error are given on the panels of this column and in Table [TJ 

To be consistent with the measured value of At, we compare our estimate of fy with 
that from the relation At/50 = (/y/200) 7 adopted in the CCF analysis. The relative 
difference of Sfv/fv = (fv — fv(At))/ fy(At) is shown in Table [TJ We see that \5fvJ fv\ 
for 7 = 0.5 is less than 10% for all the models 1 to 3, which is considerably improved over 
our previous estimate of 40% by using the simple flux mean. 

The result of fy for different values of 7 = 0.2 and 0.8 tabulated in Table [TJ shows little 
difference from that for 7 = 0.5, which indicates that estimated fy is almost independent 
of this parameter. In addition, various other light curves in asymmetric shapes also give an 
essentially same result. 
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5. Application to NGC4151 

Given that the consistent estimation of At and fy has been confirmed as above by 
using the model light curves, we apply it to the real data for NGC 4151. We divided the 
whole monitoring period into eight separate periods, each of which contains a single feature 
of flux variation, as above, with a sufficient number of more than 40 flux data pairs for 
the CCF analysis. The middle date of each period is taken as an epoch that represents 
each period. Following the procedure described in Section [31 the Monte-Carlo simulation 
was repeated 500 times for each period to form the distribution of CC centroid lag from 
which At was obtained as the median and its error as ± 34.1 percentile from the median. 
Similarly, we calculated the distribution of flux mean from which fy and its error were 
obtained. Their results for eight separate periods are given in Table [2] as well as in Figure [3j 

We here note that the A-band data of NGC 4151 were well sampled and the observed 
K light curves were smooth enough. Therefore, in order to avoid additional IR flux variation 
which consequently overestimates the errors in At, use of the SF for the interpolation was 
turned off and the simple linear interpolation was adopted for the A-band data in making 
pairs with the V^-band data. 

It is apparent that At is not constant in time. For the most extreme case, the difference 
of At between the earlier periods of 1 through 4 and the later periods of 6 through 8 is larger 
than the measurement error by many as. Since At corresponds to a light-travel distance 
to the surrounding dust from the AGN center, this should be the first co nvincing evidence 



that the inner radius of du st torus did change in an individual AGN (cf. 



19991 : 



Minezaki et al 



20041 ) 



Oknyanskij et al. 
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6. Discussion 

The change of inner radius of dust torus should occur reflecting the variation of the 
incident UV/optical flux from the central energy source. Some theoretical models suggest 



the e xistence of sublimation radius inside which dust grains are sublimated (IBarvainis 



19871 . 



19921 ). The sublimation radius expands when the UV/optical continuum becomes bright, 
then it retreats when the continuum becomes faint. If dust grains were sublimated or 
replenished immediately after the UV/optical flux variation, like the model we assumed in 
Section [31 the expected changes of At and fy should trace the simple relation of At oc L 5 
which has the slope indicated by the thick line in Figure [3j We see from this figure that 
the observed change of At follows the variation of fy in overall trend as expected. These 
results indicate that the dust torus is not a distinct, separate physical structure but is a 
part of continuous structure starting from the BLR component out to dust component with 
its inner dust edge set by the sublimation radius tha t changes with the UV/optical flux 



variation (jSuganuma et al 



2006 



Nenkova et al. 



2008|). 



On the other hand, however, the observed time-changing track does not exactly trace 
the simple relation of At oc L 5 . The lag At at the period 1 is well above that expected for 
My, and does not change so much with brightening of about 1 mag from the period 1 to 
3. Then, with subsequent fading from the period 3 to 6, At still remains constant until the 
period 4, and starts to rapidly decrease thereafter. These deviant behaviors suggest that 
the inner dust torus did not reach an equilibrium state immediately after the UV/optical 
flux variation. 

In the case of period 1, the light curves in Figure [T] indicate that there was a brighter 
period before the period 1. The expanded inner radius of dust torus at the brighter period 
could still remain larger than the radius that would be expected at the period 1. In this 
context, the behavior of almost constant At followed by its rapid decrease in the fading 
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period places a constraint such that a time scale of dust replenishment in the central region 
of AGN should be as long as the elapsed time of about one year from the period 4 through 
6. A promising mechanism of such replenishment is re-formation of dust grains rather than 
their infall from the outer region. This is understood by considering that the decrease of At 
from the period 4 to 6 is 24.5 d or a light-travel distance of 6.35 x 10 11 km in the interval 
of 309 d between the middle dates of these periods. If such decrease of At would occur 
with redistribution of dust grains supplied from the outer region of dust torus, the infall 
velocity would be 2.4 x 10 4 km s~ x or 7.9% of the light velocity. This seems highly unlikely, 
when we compare this infall velocity with a few xlO 3 km s~ x for the velocity dispersion 
of BLR clouds which exist just near the inner dust torus. Consequently, we conclude that 
re-formation of dust grains did occur in the central region where they had been sublimated. 

After the dust re-formation period, At turned almost constant again from the period 
6 through 8. This indicates that dust sublimation radius became larger than actual inner 
radius of dust torus and dust sublimation would start. The length of period of constant At 
suggests that a time scale of dust destruction could be longer than about one year. Exact 
measurement of this time scale might have been possible if the light curves were obtained 
well after the period 8. 

The formulation for improving the accuracy of lag and luminosity measurements was 
applied to the light curves for eight separate periods containing single features. This 
formulation, when applied to the light curves for the whole monitoring period containing 
several features, gives an almost average of their respective At estimates derived from single 
features, and is still useful in discussing an overall lag-luminosity relationship for a sample 
of many AGNs with a wide range of luminosity. 

We thank the staff at the Haleakala Observatories for their help with facility 
maintenance. This research has been supported partly by the Grants-in-Aid of Scientific 
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Fig. 1. — The V and K light curves of NGC 4151 (dots) in units of mJy along the vertical 
axis on the left, and the lag times At (crosses) in units of days along the vertical axis on the 
right. The horizontal bar to each cross does not specify the error but the time interval for 
which At is estimated, and the vertical bar indicates the error of At by the CCF simulations. 
The number beside each cross refers to the serial number given to each monitoring period 
in Table EJ 
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Fig. 2. — Models of V and K light curves (left column), and the corresponding distributions 
of CC centroid lag (middle column) and mean V flux (right column), for the cases of time- 
symmetric incident V light curve. Three panels from left to right in each row shows a 
procedure of estimating At and fy from each model of light curve (see text). The solid lines 
on the panels of middle column show the CCFs with their values on the right axis. The 
horizontal arrows on the panels of left column indicate the periods for which the flux data 
are excluded for estimation of fy. Note that the units are arbitrary in flux. 
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Fig. 3. — The time-changing track of NGC 4151 in the lag time versus absolute ^-magnitude 
diagram, as shown by the dot-connecting thin line. The number beside each dot refers to 
the serial number given to each monitoring period in Table 2. The thick line is the scaling 
relation of At oc L 5 expected from dust reverberation. We used H = 70km s -1 Mpc -1 for 
estimation of My. 
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Table 1. Results of simulation runs for model light curves 



Model 


Type 


'y 


At ( days) 


fv 


Sfv/fv 


1 
2 
3 


Linear bottom 
Quadratic bottom 
Gaussian peak 


0.5 
0.5 
0.5 


74.6li; 6 7 
80.3±1;S 

' '-°-0.8 


483.6±|"i 
482.9±f , j 
488.911;-* 


0.09 
-0.06 
0.02 


la 
2a 


Linear bottom 
Quadratic bottom 


0.2 
0.2 


59.91JJ 


489.4lfj 
475.21';; 


0.06 
-0.04 


lb 
2b 


Linear bottom 
Quadratic bottom 


0.8 
0.8 


107.3^'j 
107.1±^;| 


490.91^ 
481.6±™ 


-0.06 
-0.05 



Note. — 7 is a power index for the generator At/50 = (/V/200) 7 . Sfv/fv 
is defined as (fv — /y(Af))//y (At), where fv(At) is estimated from the 
generator with the measured value of At. The units here are days for At and 
arbitrary for fv- 



Table 2. Lag time and ^-magnitude of NGC 4151 for separate monitoring periods 



Period MJD (days) At (days) my (mag) 



1 


51915 


.6-52064. 


3 


58. 


,4 


+1.1 

-1.2 


14. 


,49 


+0.05 
-0.06 


2 


52265 


.5-52643. 


6 


71, 


1 


+ 13.5 
-9.8 


13. 


,04 


+0.07 
-0.07 


3 


52610. 


.6-52967. 


6 


60 


.6 


+ 1.5 
-1.8 


12 


.69 


+0.04 

-0.04 


4 


53039 


.4-53356 


5 


65 


1 


+2.9 
-3.0 


13 


,50 


+0.07 
-0.08 


5 


53148 


.4-53490 


3 


50. 


.4 


+6.9 
-12.1 


13 


,66 


+0.09 
-0.10 


6 


53436. 


.4-53578 


3 


40, 


6 


+ 1.5 
-1.3 


14 


15 


+0.05 

-0.06 


7 


53687. 


,6-53864. 


.4 


33 


1 


+2.3 
-2.3 


13 


,54 


+0.06 
-0.06 


8 


53753. 


,4-53921. 


3 


38, 


.5 


+ 1.3 
-1.2 


13 


,70 


+0.10 
-0.10 



